library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ stringr 1.4.0
## ✓ tidyr   1.1.4     ✓ forcats 0.5.1
## ✓ readr   2.1.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(knitr)
library(ggplot2)
library(skimr)
f<- "https://raw.githubusercontent.com/mrpickett26/Helping-Maddie/main/M%20Divine%20VC%20Psychfest%20Data.csv"
d<- read_csv(f, col_names=TRUE)
## New names:
## * `Participant ID` -> `Participant ID...1`
## * `Start Date` -> `Start Date...56`
## * `End Date` -> `End Date...57`
## * Progress -> Progress...58
## * `Duration (in seconds)` -> `Duration (in seconds)...59`
## * ...
## Rows: 210 Columns: 256
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (26): Semester, Hormones D1 Complete, Hormones D2 Complete, VAS D1 Comp...
## dbl (227): Participant ID...1, Sex, BirthControl, Do you take any chemical c...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d <-d %>% dplyr::select(Sex,BirthControl,D2S1_Cortisol,D2S1_Testosterone,D2S1_Progesterone,D2S1_Estradiol,D2S2_Cortisol,D2S2_Testosterone,D2S2_Progesterone,D2S2_Estradiol,D2S3_Cortisol,D2S3_Testosterone,D2S3_Progesterone,D2S3_Estradiol,D2VAS1_Stress,D2VAS1_Shame,D2VAS2_Stress,D2VAS2_Shame,D2VAS3_Stress,D2VAS3_Shame,D2VAS4_Stress,D2VAS4_Shame,D2VAS5_Stress,D2VAS5_Shame,CESD)
d %>% group_by(Sex) %>% summarise_all(~min(.x,na.rm=TRUE))
## Warning in min(.x, na.rm = TRUE): no non-missing arguments to min; returning Inf
## # A tibble: 2 × 25
##     Sex BirthControl D2S1_Cortisol D2S1_Testosterone D2S1_Progesterone
##   <dbl>        <dbl>         <dbl>             <dbl>             <dbl>
## 1     1          Inf         0.103              8.02             0.470
## 2     2            0         0.353              1.29             0.769
## # … with 20 more variables: D2S1_Estradiol <dbl>, D2S2_Cortisol <dbl>,
## #   D2S2_Testosterone <dbl>, D2S2_Progesterone <dbl>, D2S2_Estradiol <dbl>,
## #   D2S3_Cortisol <dbl>, D2S3_Testosterone <dbl>, D2S3_Progesterone <dbl>,
## #   D2S3_Estradiol <dbl>, D2VAS1_Stress <dbl>, D2VAS1_Shame <dbl>,
## #   D2VAS2_Stress <dbl>, D2VAS2_Shame <dbl>, D2VAS3_Stress <dbl>,
## #   D2VAS3_Shame <dbl>, D2VAS4_Stress <dbl>, D2VAS4_Shame <dbl>,
## #   D2VAS5_Stress <dbl>, D2VAS5_Shame <dbl>, CESD <dbl>
d<- d %>% mutate(Sex = Sex -1)
d<-d %>% mutate(BirthControl = replace_na(BirthControl,0))
data_v1<- d %>% mutate(
  #adjusts minimum for each gender to be 1 by shifting all values over
  across(ends_with("Cortisol"),~ifelse(Sex == 1, .+(1-0.1031219), .+(1-0.3446964))),
  across(ends_with("Progesterone"),~ifelse(Sex == 1, .+(1- 0.46972), .+(1-0.68095))),
  across(ends_with("Estradiol"),~ifelse(Sex == 1, .+(1- 0.31000  ), .+(1-0.31622)))
  ) %>% mutate(
    #takes the log of all the hormones for all 3 periods in D2
    across(starts_with("D2S"),~log(.))
  )

fn_auc_g = function(s1,s2,s3){
  (s1+s2)*40/2 + (s2+s3)*15/2
}
fn_auc_i = function(s1,s2,s3){
  (s1+s2)*40/2 + (s2+s3)*15/2 - s1*(55)
}

data_v1 <-data_v1 %>% mutate(D2_Cortisol_AUC_i = fn_auc_i(D2S1_Cortisol,D2S2_Cortisol,D2S3_Cortisol),
                   D2_Cortisol_AUC_g = fn_auc_g(D2S1_Cortisol,D2S2_Cortisol,D2S3_Cortisol),
                   D2_Testosterone_AUC_i = fn_auc_i(D2S1_Testosterone,D2S2_Testosterone,D2S3_Testosterone),
                   D2_Testosterone_AUC_g = fn_auc_g(D2S1_Testosterone,D2S2_Testosterone,D2S3_Testosterone),
                   D2_Progesterone_AUC_i = fn_auc_i(D2S1_Progesterone,D2S2_Progesterone,D2S3_Progesterone),
                   D2_Progesterone_AUC_g = fn_auc_g(D2S1_Progesterone,D2S2_Progesterone,D2S3_Progesterone),
                   D2_Estradiol_AUC_i = fn_auc_i(D2S1_Estradiol,D2S2_Estradiol,D2S3_Estradiol),
                   D2_Estradiol_AUC_g = fn_auc_g(D2S1_Estradiol,D2S2_Estradiol,D2S3_Estradiol)
                   )

hist(data_v1$CESD)

hist(data_v1$D2_Cortisol_AUC_i)

hist(data_v1$D2_Cortisol_AUC_g)

hist(data_v1$D2_Progesterone_AUC_i)

hist(data_v1$D2_Progesterone_AUC_g)

hist(data_v1$D2_Testosterone_AUC_i)

hist(data_v1$D2_Estradiol_AUC_i)

hist(data_v1$D2_Estradiol_AUC_g)

# library(GGally)
# ggpairs(data_v1 %>% filter(Sex==1), columns = c("CESD","D2_Cortisol_AUC_i","D2_Testosterone_AUC_i","D2_Estradiol_AUC_i","D2_Progesterone_AUC_i"))
# 
# ggpairs(data_v1 %>% filter(Sex==0), columns = c("CESD","D2_Cortisol_AUC_i","D2_Testosterone_AUC_i","D2_Estradiol_AUC_i","D2_Progesterone_AUC_i"))
# 
# ggpairs(data_v1 %>% filter(Sex==1), columns = c("CESD","D2_Cortisol_AUC_g","D2_Testosterone_AUC_g","D2_Estradiol_AUC_g","D2_Progesterone_AUC_g"))

# Cesd predicted by sex
m1<- lm(data=data_v1, CESD~Sex)
summary(m1)
## 
## Call:
## lm(formula = CESD ~ Sex, data = data_v1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.076  -7.076  -2.220   5.924  28.636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  35.0764     0.7937   44.19   <2e-16 ***
## Sex          -1.7128     1.4158   -1.21    0.228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.525 on 208 degrees of freedom
## Multiple R-squared:  0.006987,   Adjusted R-squared:  0.002213 
## F-statistic: 1.463 on 1 and 208 DF,  p-value: 0.2278
data_female<- data_v1%>%filter(Sex==1)
data_male<- data_v1%>% filter(Sex==0)
# p value greater than 0.05 so no significant difference 
relCESD<- lm(CESD~BirthControl, data=data_female, na.action = na.exclude)

## Now going to run the full linear model for females, accounting for the interaction effect of  birth control for all predictiors, as well as the anova on the linear model 
library(broom)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
## 
##     some
## The following object is masked from 'package:dplyr':
## 
##     recode
lm_cortisol_bc_i<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_i+D2_Testosterone_AUC_i+D2_Progesterone_AUC_i+D2_Estradiol_AUC_i)))
summary(lm_cortisol_bc_i)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_i + D2_Testosterone_AUC_i + 
##     D2_Progesterone_AUC_i + D2_Estradiol_AUC_i)), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.755  -4.479  -1.179   3.730  22.530 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        32.78747    1.43076  22.916  < 2e-16 ***
## BirthControl                        1.39401    2.69321   0.518  0.60702    
## D2_Cortisol_AUC_i                   0.06606    0.05693   1.160  0.25136    
## D2_Testosterone_AUC_i              -0.02858    0.06972  -0.410  0.68360    
## D2_Progesterone_AUC_i              -0.13369    0.12671  -1.055  0.29646    
## D2_Estradiol_AUC_i                  0.16577    0.11068   1.498  0.14050    
## BirthControl:D2_Cortisol_AUC_i      0.05258    0.10600   0.496  0.62208    
## BirthControl:D2_Testosterone_AUC_i  0.10690    0.11422   0.936  0.35382    
## BirthControl:D2_Progesterone_AUC_i  0.41069    0.16875   2.434  0.01856 *  
## BirthControl:D2_Estradiol_AUC_i    -0.66104    0.20976  -3.151  0.00274 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.053 on 50 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.3418, Adjusted R-squared:  0.2234 
## F-statistic: 2.885 on 9 and 50 DF,  p-value: 0.007966
plot(fitted(lm_cortisol_bc_i), residuals(lm_cortisol_bc_i))

hist(residuals(lm_cortisol_bc_i))

qqnorm(residuals(lm_cortisol_bc_i))
qqline(residuals(lm_cortisol_bc_i))

m.aov_i <- Anova(lm_cortisol_bc_i, type = "II")
m.aov_i
## Anova Table (Type II tests)
## 
## Response: CESD
##                                    Sum Sq Df F value   Pr(>F)   
## BirthControl                        120.5  1  1.8580 0.178968   
## D2_Cortisol_AUC_i                   185.5  1  2.8610 0.096974 . 
## D2_Testosterone_AUC_i                 2.7  1  0.0415 0.839424   
## D2_Progesterone_AUC_i                88.7  1  1.3672 0.247841   
## D2_Estradiol_AUC_i                    2.5  1  0.0378 0.846575   
## BirthControl:D2_Cortisol_AUC_i       16.0  1  0.2460 0.622079   
## BirthControl:D2_Testosterone_AUC_i   56.8  1  0.8759 0.353817   
## BirthControl:D2_Progesterone_AUC_i  384.1  1  5.9229 0.018561 * 
## BirthControl:D2_Estradiol_AUC_i     644.1  1  9.9317 0.002744 **
## Residuals                          3242.7 50                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## In this case we see significance in differences in CESD scores from cortisol as a predictor alone (0.9 sig), Progesterone with BC accounting for interaction effects (0.95), and Estradiol accounting for BC interaction effects (0.99)

#Now I will do the same for the AUCg, sorry the variables are fucked 
lm_cortisol_bc_g<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g)))
summary(lm_cortisol_bc_g)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_g + D2_Testosterone_AUC_g + 
##     D2_Progesterone_AUC_g + D2_Estradiol_AUC_g)), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.470  -5.701  -1.551   3.481  23.033 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         41.25707   10.66338   3.869 0.000317 ***
## BirthControl                       -33.11539   15.30073  -2.164 0.035243 *  
## D2_Cortisol_AUC_g                   -0.05453    0.04594  -1.187 0.240827    
## D2_Testosterone_AUC_g                0.06216    0.05942   1.046 0.300568    
## D2_Progesterone_AUC_g               -0.02303    0.01732  -1.329 0.189739    
## D2_Estradiol_AUC_g                  -0.07103    0.11128  -0.638 0.526162    
## BirthControl:D2_Cortisol_AUC_g       0.07432    0.05491   1.354 0.181964    
## BirthControl:D2_Testosterone_AUC_g  -0.02272    0.08946  -0.254 0.800563    
## BirthControl:D2_Progesterone_AUC_g   0.09045    0.04902   1.845 0.070957 .  
## BirthControl:D2_Estradiol_AUC_g      0.33400    0.16912   1.975 0.053818 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.458 on 50 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.274,  Adjusted R-squared:  0.1433 
## F-statistic: 2.097 on 9 and 50 DF,  p-value: 0.04746
plot(fitted(lm_cortisol_bc_g), residuals(lm_cortisol_bc_g))

hist(residuals(lm_cortisol_bc_g))

qqnorm(residuals(lm_cortisol_bc_g))
qqline(residuals(lm_cortisol_bc_g))

m.aov_g <- Anova(lm_cortisol_bc_g, type = "II")
m.aov_g
## Anova Table (Type II tests)
## 
## Response: CESD
##                                    Sum Sq Df F value  Pr(>F)  
## BirthControl                         52.1  1  0.7288 0.39733  
## D2_Cortisol_AUC_g                     0.7  1  0.0099 0.92097  
## D2_Testosterone_AUC_g                98.6  1  1.3776 0.24607  
## D2_Progesterone_AUC_g                37.5  1  0.5243 0.47238  
## D2_Estradiol_AUC_g                   55.1  1  0.7706 0.38422  
## BirthControl:D2_Cortisol_AUC_g      131.1  1  1.8321 0.18196  
## BirthControl:D2_Testosterone_AUC_g    4.6  1  0.0645 0.80056  
## BirthControl:D2_Progesterone_AUC_g  243.5  1  3.4042 0.07096 .
## BirthControl:D2_Estradiol_AUC_g     279.0  1  3.9001 0.05382 .
## Residuals                          3576.9 50                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##In this case we see significance in differences in Estradiol with BC for CESD score response accounting for interaction effects (0.9), and Progesterone accounting for BC interaction effects (0.99)

##Now moving on to Male AUCg
lm_male_g<-lm(data=data_male,CESD~(D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g))
summary(lm_male_g)
## 
## Call:
## lm(formula = CESD ~ (D2_Cortisol_AUC_g + D2_Testosterone_AUC_g + 
##     D2_Progesterone_AUC_g + D2_Estradiol_AUC_g), data = data_male)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.679  -6.884  -1.992   5.218  23.004 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           34.639800   8.466501   4.091 7.95e-05 ***
## D2_Cortisol_AUC_g     -0.017297   0.034339  -0.504    0.615    
## D2_Testosterone_AUC_g  0.015481   0.036185   0.428    0.670    
## D2_Progesterone_AUC_g -0.004648   0.028829  -0.161    0.872    
## D2_Estradiol_AUC_g    -0.011616   0.071309  -0.163    0.871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.633 on 116 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.003139,   Adjusted R-squared:  -0.03124 
## F-statistic: 0.09131 on 4 and 116 DF,  p-value: 0.985
plot(fitted(lm_male_g), residuals(lm_male_g))

hist(residuals(lm_male_g))

qqnorm(residuals(lm_male_g))
qqline(residuals(lm_male_g))

male.aov_g <- Anova(lm_male_g, type = "II")
male.aov_g
## Anova Table (Type II tests)
## 
## Response: CESD
##                        Sum Sq  Df F value Pr(>F)
## D2_Cortisol_AUC_g        23.5   1  0.2537 0.6154
## D2_Testosterone_AUC_g    17.0   1  0.1830 0.6696
## D2_Progesterone_AUC_g     2.4   1  0.0260 0.8722
## D2_Estradiol_AUC_g        2.5   1  0.0265 0.8709
## Residuals             10765.1 116
##No reported significant difference in Hormone Levels with respect to CESD scores

#Now moving on to Malue AUCi
lm_male_i<-lm(data=data_male,CESD~(D2_Cortisol_AUC_i+D2_Testosterone_AUC_i+D2_Progesterone_AUC_i+D2_Estradiol_AUC_i))
summary(lm_male_i)
## 
## Call:
## lm(formula = CESD ~ (D2_Cortisol_AUC_i + D2_Testosterone_AUC_i + 
##     D2_Progesterone_AUC_i + D2_Estradiol_AUC_i), data = data_male)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.753  -6.793  -2.485   5.806  23.303 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           35.28293    1.10938  31.804   <2e-16 ***
## D2_Cortisol_AUC_i      0.01094    0.03752   0.292    0.771    
## D2_Testosterone_AUC_i  0.03135    0.08092   0.387    0.699    
## D2_Progesterone_AUC_i -0.02998    0.05809  -0.516    0.607    
## D2_Estradiol_AUC_i    -0.06328    0.10349  -0.611    0.542    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.614 on 116 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.007213,   Adjusted R-squared:  -0.02702 
## F-statistic: 0.2107 on 4 and 116 DF,  p-value: 0.932
plot(fitted(lm_male_i), residuals(lm_male_i))

hist(residuals(lm_male_i))

qqnorm(residuals(lm_male_i))
qqline(residuals(lm_male_i))

male.aov_i <- Anova(lm_male_i, type = "II")
male.aov_i
## Anova Table (Type II tests)
## 
## Response: CESD
##                        Sum Sq  Df F value Pr(>F)
## D2_Cortisol_AUC_i         7.9   1  0.0850 0.7711
## D2_Testosterone_AUC_i    13.9   1  0.1501 0.6992
## D2_Progesterone_AUC_i    24.6   1  0.2664 0.6068
## D2_Estradiol_AUC_i       34.6   1  0.3738 0.5421
## Residuals             10721.1 116
##No reported significant difference in Hormone Levels with respect to CESD scores


##Look at stressors with respect to hormones 
lm_male_i_stress<-lm(data=data_male,CESD~(D2VAS1_Stress+D2VAS2_Stress+D2VAS3_Stress+D2VAS4_Stress+D2VAS5_Stress))
summary(lm_male_i_stress)
## 
## Call:
## lm(formula = CESD ~ (D2VAS1_Stress + D2VAS2_Stress + D2VAS3_Stress + 
##     D2VAS4_Stress + D2VAS5_Stress), data = data_male)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.683  -6.548  -1.816   6.273  21.871 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.53111    1.50543  18.952   <2e-16 ***
## D2VAS1_Stress  0.07294    0.06222   1.172   0.2433    
## D2VAS2_Stress  0.08730    0.04349   2.007   0.0468 *  
## D2VAS3_Stress  0.03547    0.04704   0.754   0.4522    
## D2VAS4_Stress  0.08163    0.11033   0.740   0.4608    
## D2VAS5_Stress -0.05930    0.11953  -0.496   0.6207    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.938 on 127 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.1863, Adjusted R-squared:  0.1542 
## F-statistic: 5.814 on 5 and 127 DF,  p-value: 7.208e-05
plot(fitted(lm_male_i_stress), residuals(lm_male_i_stress))

hist(residuals(lm_male_i_stress))

qqnorm(residuals(lm_male_i_stress))
qqline(residuals(lm_male_i_stress))

male.aov_i_stress <- Anova(lm_male_i_stress, type = "II")
male.aov_i_stress
## Anova Table (Type II tests)
## 
## Response: CESD
##                Sum Sq  Df F value  Pr(>F)  
## D2VAS1_Stress   109.8   1  1.3743 0.24327  
## D2VAS2_Stress   321.9   1  4.0297 0.04683 *
## D2VAS3_Stress    45.4   1  0.5687 0.45218  
## D2VAS4_Stress    43.7   1  0.5474 0.46076  
## D2VAS5_Stress    19.7   1  0.2461 0.62069  
## Residuals     10145.6 127                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#So for males only timepoint 2 is significant for differences in CESD scores

##Now do this for Shame 
lm_male_i_shame<-lm(data=data_male,CESD~(D2VAS1_Shame+D2VAS2_Shame+D2VAS3_Shame+D2VAS4_Shame+D2VAS5_Shame))
summary(lm_male_i_shame)
## 
## Call:
## lm(formula = CESD ~ (D2VAS1_Shame + D2VAS2_Shame + D2VAS3_Shame + 
##     D2VAS4_Shame + D2VAS5_Shame), data = data_male)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.888  -7.430  -2.008   6.538  22.252 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  32.37114    1.03619  31.240   <2e-16 ***
## D2VAS1_Shame  0.06015    0.10131   0.594   0.5538    
## D2VAS2_Shame -0.06769    0.06355  -1.065   0.2888    
## D2VAS3_Shame  0.11173    0.04996   2.236   0.0271 *  
## D2VAS4_Shame  0.09076    0.10034   0.904   0.3674    
## D2VAS5_Shame -0.07813    0.10906  -0.716   0.4750    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.366 on 127 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.1064, Adjusted R-squared:  0.07121 
## F-statistic: 3.024 on 5 and 127 DF,  p-value: 0.01298
plot(fitted(lm_male_i_shame), residuals(lm_male_i_shame))

hist(residuals(lm_male_i_shame))

qqnorm(residuals(lm_male_i_shame))
qqline(residuals(lm_male_i_shame))

male.aov_i_shame <- Anova(lm_male_i_shame, type = "II")
male.aov_i_shame
## Anova Table (Type II tests)
## 
## Response: CESD
##               Sum Sq  Df F value  Pr(>F)  
## D2VAS1_Shame    30.9   1  0.3525 0.55376  
## D2VAS2_Shame    99.5   1  1.1345 0.28883  
## D2VAS3_Shame   438.7   1  5.0007 0.02708 *
## D2VAS4_Shame    71.8   1  0.8181 0.36745  
## D2VAS5_Shame    45.0   1  0.5133 0.47505  
## Residuals    11141.3 127                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#So for males only timepoint 3 is significant for shame for differences in CESD scores

##Do the same for females
lm_female_i_stress<-lm(data=data_female,CESD~(D2VAS1_Stress+D2VAS2_Stress+D2VAS3_Stress+D2VAS4_Stress+D2VAS5_Stress))
summary(lm_female_i_stress)
## 
## Call:
## lm(formula = CESD ~ (D2VAS1_Stress + D2VAS2_Stress + D2VAS3_Stress + 
##     D2VAS4_Stress + D2VAS5_Stress), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.821  -4.527  -1.060   4.895  20.484 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.2051875  2.0130020  14.012   <2e-16 ***
## D2VAS1_Stress -0.0002214  0.0710737  -0.003   0.9975    
## D2VAS2_Stress  0.0011820  0.0564798   0.021   0.9834    
## D2VAS3_Stress  0.1284010  0.0524936   2.446   0.0181 *  
## D2VAS4_Stress -0.0509880  0.0944010  -0.540   0.5916    
## D2VAS5_Stress -0.0549777  0.0966732  -0.569   0.5722    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.192 on 49 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.1824, Adjusted R-squared:  0.09896 
## F-statistic: 2.186 on 5 and 49 DF,  p-value: 0.07081
plot(fitted(lm_female_i_stress), residuals(lm_female_i_stress))

hist(residuals(lm_female_i_stress))

qqnorm(residuals(lm_female_i_stress))
qqline(residuals(lm_female_i_stress))

female.aov_i_stress <- Anova(lm_female_i_stress, type = "II")
female.aov_i_stress
## Anova Table (Type II tests)
## 
## Response: CESD
##                Sum Sq Df F value  Pr(>F)  
## D2VAS1_Stress    0.00  1  0.0000 0.99753  
## D2VAS2_Stress    0.02  1  0.0004 0.98339  
## D2VAS3_Stress  309.45  1  5.9831 0.01808 *
## D2VAS4_Stress   15.09  1  0.2917 0.59156  
## D2VAS5_Stress   16.73  1  0.3234 0.57216  
## Residuals     2534.35 49                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#So for females only timepoint 3 (stress) is significant for differences in CESD scores

#shame for females
lm_female_i_shame<-lm(data=data_female,CESD~(D2VAS1_Shame+D2VAS2_Shame+D2VAS3_Shame+D2VAS4_Shame+D2VAS5_Shame))
summary(lm_female_i_shame)
## 
## Call:
## lm(formula = CESD ~ (D2VAS1_Shame + D2VAS2_Shame + D2VAS3_Shame + 
##     D2VAS4_Shame + D2VAS5_Shame), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.794  -4.596  -1.396   3.267  22.467 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.6598992  1.3169998  21.762  < 2e-16 ***
## D2VAS1_Shame  0.0735331  0.1513691   0.486 0.629284    
## D2VAS2_Shame  0.0002647  0.0742381   0.004 0.997170    
## D2VAS3_Shame  0.1479441  0.0406896   3.636 0.000664 ***
## D2VAS4_Shame -0.1725549  0.1129552  -1.528 0.133032    
## D2VAS5_Shame  0.0367721  0.1270935   0.289 0.773549    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.869 on 49 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.2541, Adjusted R-squared:  0.1779 
## F-statistic: 3.338 on 5 and 49 DF,  p-value: 0.01132
plot(fitted(lm_female_i_shame), residuals(lm_female_i_shame))

hist(residuals(lm_female_i_shame))

qqnorm(residuals(lm_female_i_shame))
qqline(residuals(lm_female_i_shame))

female.aov_i_shame <- Anova(lm_female_i_shame, type = "II")
female.aov_i_shame
## Anova Table (Type II tests)
## 
## Response: CESD
##               Sum Sq Df F value   Pr(>F)    
## D2VAS1_Shame   11.14  1  0.2360 0.629284    
## D2VAS2_Shame    0.00  1  0.0000 0.997170    
## D2VAS3_Shame  623.82  1 13.2199 0.000664 ***
## D2VAS4_Shame  110.12  1  2.3337 0.133032    
## D2VAS5_Shame    3.95  1  0.0837 0.773549    
## Residuals    2312.21 49                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Some things I tried out... didnt like it as much
# Male_Estradiol_g<- data_male$D2_Estradiol_AUC_g
# Male_Test_g<- data_male$D2_Testosterone_AUC_g
# Male_Cort_g<- data_male$D2_Cortisol_AUC_g
# Male_Proges_g<- data_male$D2_Progesterone_AUC_g
# Male_Estradiol_i<- data_male$D2_Estradiol_AUC_i
# Male_Test_i<- data_male$D2_Testosterone_AUC_i
# Male_Cort_i<- data_male$D2_Cortisol_AUC_i
# Male_Proges_i<- data_male$D2_Progesterone_AUC_i
# Male_Estradiol_g_norm <- rnorm(200,mean=mean(Male_Estradiol_g, na.rm=TRUE), sd=sd(Male_Estradiol_g, na.rm=TRUE))
# Male_Test_g_norm <- rnorm(200,mean=mean(Male_Test_g, na.rm=TRUE), sd=sd(Male_Test_g, na.rm=TRUE))
# Male_Cort_g_norm <- rnorm(200,mean=mean(Male_Cort_g, na.rm=TRUE), sd=sd(Male_Cort_g, na.rm=TRUE))
# Male_Proges_g_norm <- rnorm(200,mean=mean(Male_Proges_g, na.rm=TRUE), sd=sd(Male_Proges_g, na.rm=TRUE))
# Male_Estradiol_i_norm <- rnorm(200,mean=mean(Male_Estradiol_i, na.rm=TRUE), sd=sd(Male_Estradiol_g, na.rm=TRUE))
# Male_Test_i_norm <- rnorm(200,mean=mean(Male_Test_i, na.rm=TRUE), sd=sd(Male_Test_g, na.rm=TRUE))
# Male_Cort_i_norm <- rnorm(200,mean=mean(Male_Cort_i, na.rm=TRUE), sd=sd(Male_Cort_g, na.rm=TRUE))
# Male_Proges_i_norm <- rnorm(200,mean=mean(Male_Proges_i, na.rm=TRUE), sd=sd(Male_Proges_g, na.rm=TRUE))
# 
# plot(x=data_male$CESD, y=Male_Estradiol_g,Male_Estradiol_i, Male_Test_g,Male_Test_i, Male_Cort_g, Male_Cort_i,Male_Proges_g,Male_Proges_i, na.rm=FALSE)
# main = "Male Hormone Levels",
# at = c(1,2,4,5),
# names = c("AUCg Estradiol", "AUCi Estradiol", "AUCg Testosterone","AUCi Testosterone", "AUCg Cortisol","AUCi Cortisol", "AUCg Progresterone","AUCi Progresterone"),
# las = 2,
# col = c("orange","red"),
# border = "brown",
# horizontal = FALSE,
# notch = TRUE


# data_female<-data_female%>%mutate(
#   resid_CESD=residuals(relCESD)
# )
# 
# m_null<-lm(data=data_female, resid_CESD~1)
# summary(m_null)
# 
# add1(m_null, scope=.~.+D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g, test="F")
# 
# m2<-update(m_null, formula=.~. +D2_Testosterone_AUC_g)
# add1(m2, scope=.~.+D2_Cortisol_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g, test="F")
# 
# m3<-update(m2, formula=.~. +D2_Estradiol_AUC_g)
# add1(m3, scope=.~.+D2_Cortisol_AUC_g+D2_Progesterone_AUC_g, test="F")
# m4<-update(m3, formula=.~. +D2_Progesterone_AUC_g)
# add1(m4, scope=.~.+D2_Cortisol_AUC_g, test="F")
# 
# 
#   lm(data=data_female, resid_CESD~D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g)
# summary(m_full)
# 
# 
# m2<- lm(data=data_female, CESD~BirthControl)
# summary(m2)
# 
# m3<- lm(data=data_female, CESD~BirthControl+)
# summary(m3)
# 
# m2<- lm(data=data_v1, CESD~ logRange+Migration)
# m3<- lm(data=data_v1, CESD~logRange)
# m4<- lm(data=data_v1, CESD~Migration)
# m5<- lm(data=d, CESD~1)
# 
# anova(m1,m2, test="F")

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

#So for females only timepoint 3 (shame) is significant for differences in CESD scores, VERY signifcant (want to emphasize this 0.999)
library(sjPlot)
## Registered S3 method overwritten by 'parameters':
##   method                         from      
##   format.parameters_distribution datawizard
library(sjmisc)
## 
## Attaching package: 'sjmisc'
## The following object is masked from 'package:skimr':
## 
##     to_long
## The following object is masked from 'package:purrr':
## 
##     is_empty
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## The following object is masked from 'package:tibble':
## 
##     add_case
library(ggplot2)
theme_set(theme_sjplot())


#Now it only makes sense to look at hormonal female data and discuss it, run the signigicant linear models

#Do this for the female AUCi significant data
lm_cortisol_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_i)))
summary(lm_cortisol_i_sig)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_i)), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.666  -5.280  -1.236   5.225  24.235 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    32.78410    1.49281  21.961   <2e-16 ***
## BirthControl                    2.03622    2.32926   0.874   0.3857    
## D2_Cortisol_AUC_i               0.04766    0.05485   0.869   0.3886    
## BirthControl:D2_Cortisol_AUC_i  0.20009    0.09646   2.074   0.0427 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.562 on 56 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.1668, Adjusted R-squared:  0.1222 
## F-statistic: 3.737 on 3 and 56 DF,  p-value: 0.01611
plot(lm_cortisol_i_sig)

plot_model(lm_cortisol_i_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_cortisol_i_sig, type = "std")

lm_estradiol_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Estradiol_AUC_i)))
summary(lm_estradiol_i_sig)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Estradiol_AUC_i)), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.001  -6.034  -2.595   4.894  23.816 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      33.0233     1.4848  22.240   <2e-16 ***
## BirthControl                      0.8283     2.5188   0.329   0.7435    
## D2_Estradiol_AUC_i                0.1575     0.1195   1.319   0.1926    
## BirthControl:D2_Estradiol_AUC_i  -0.4383     0.2110  -2.077   0.0424 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.988 on 56 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.0818, Adjusted R-squared:  0.03261 
## F-statistic: 1.663 on 3 and 56 DF,  p-value: 0.1854
plot(lm_estradiol_i_sig)

plot_model(lm_estradiol_i_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_estradiol_i_sig, type = "std")

lm_proges_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Progesterone_AUC_i)))
summary(lm_proges_i_sig)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Progesterone_AUC_i)), 
##     data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.506  -6.803  -1.996   5.351  23.929 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        33.34203    1.46483  22.762   <2e-16 ***
## BirthControl                        3.90971    2.49791   1.565   0.1232    
## D2_Progesterone_AUC_i              -0.07227    0.12008  -0.602   0.5497    
## BirthControl:D2_Progesterone_AUC_i  0.31309    0.15362   2.038   0.0463 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.819 on 56 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.116,  Adjusted R-squared:  0.06863 
## F-statistic: 2.449 on 3 and 56 DF,  p-value: 0.07303
plot(lm_proges_i_sig)

plot_model(lm_proges_i_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_proges_i_sig, type = "std")

#Do Do this for the female AUCg significant data

lm_estradiol_g_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Estradiol_AUC_g)))
summary(lm_estradiol_g_sig)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Estradiol_AUC_g)), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.413  -6.905  -1.247   4.134  24.811 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      42.4875     8.9621   4.741  1.5e-05 ***
## BirthControl                    -27.0616    12.2387  -2.211   0.0311 *  
## D2_Estradiol_AUC_g               -0.1178     0.1124  -1.048   0.2990    
## BirthControl:D2_Estradiol_AUC_g   0.4133     0.1658   2.493   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.798 on 56 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.1201, Adjusted R-squared:  0.07297 
## F-statistic: 2.548 on 3 and 56 DF,  p-value: 0.06497
plot(lm_estradiol_g_sig)

plot_model(lm_estradiol_g_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_estradiol_g_sig, type = "std")

lm_proges_g_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Progesterone_AUC_g)))
summary(lm_proges_g_sig)
## 
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Progesterone_AUC_g)), 
##     data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.005  -5.933  -2.128   5.392  23.547 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        37.10804    3.10247  11.961  < 2e-16 ***
## BirthControl                       -8.23628    4.18881  -1.966  0.05423 .  
## D2_Progesterone_AUC_g              -0.02403    0.01706  -1.408  0.16458    
## BirthControl:D2_Progesterone_AUC_g  0.12209    0.03809   3.206  0.00223 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.576 on 56 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.164,  Adjusted R-squared:  0.1192 
## F-statistic: 3.661 on 3 and 56 DF,  p-value: 0.0176
plot(lm_proges_g_sig)

plot_model(lm_proges_g_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_proges_g_sig, type = "std")

## Interesting to look at interaction effect of shame and stress of females on CESD at timepoint 3, run a linear model to do that 
lm_tress_shame_sig<-lm(data=data_female,CESD~(D2VAS3_Shame*D2VAS3_Stress))
summary(lm_tress_shame_sig)
## 
## Call:
## lm(formula = CESD ~ (D2VAS3_Shame * D2VAS3_Stress), data = data_female)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.110  -3.855  -1.314   3.252  23.450 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                30.3822801  1.8244787  16.653  < 2e-16 ***
## D2VAS3_Shame               -0.0548903  0.0535437  -1.025  0.31013    
## D2VAS3_Stress              -0.0450457  0.0496274  -0.908  0.36832    
## D2VAS3_Shame:D2VAS3_Stress  0.0027281  0.0009827   2.776  0.00768 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.528 on 51 degrees of freedom
##   (11 observations deleted due to missingness)
## Multiple R-squared:  0.2988, Adjusted R-squared:  0.2576 
## F-statistic: 7.245 on 3 and 51 DF,  p-value: 0.0003857
plot(lm_tress_shame_sig)

plot_model(lm_tress_shame_sig, type = "diag")
## [[1]]

## 
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

plot_model(lm_tress_shame_sig, type = "std")

##Interaction is significant (0.05)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.